# This file is part of Sonedyan.
#
# Sonedyan is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public
# License as published by the Free Software Foundation;
# either version 3 of the License, or (at your option) any
# later version.
#
# Sonedyan is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied
# warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
# PURPOSE.  See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public.
# If not, see <http://www.gnu.org/licenses/>.
#
# Copyright (C) 2009-2012 Jimmy Dubuisson <jimmy.dubuisson@gmail.com>

# needed for operator %^%
library(expm)

# add self loops to the adjancency matrix
addSelfLoops <- function(m)
{
  m[col(m) == row(m)] <- 1
  return(m)
}

# normalize the matrix columns
normalizeColumns <- function(m)
{
  # return(sweep(m, 2, colSums(m), "/"))
  return(scale(m, center=FALSE, scale=colSums(m)))
}

# add self loops to the adjacency matrix
# and normalize m columns
init <- function(m)
{
  m <- addSelfLoops(m)
  m <- normalizeColumns(m)
  return(m)
}

# compute col(m)^r and normalize the new matrix columns
inflation <- function(m, r)
{
  return(normalizeColumns(sweep(m, 2, r, "^")))
}

# compute m^e
expansion <- function(m, e)
{
  return(m %^% e)
}

# MCL algorithm
mcl <- function(m, e, r, n)
{
  print("Init")	
  print(m)

  m <- init(m)
  
  print("Ready")
  print(m)
 
  for (i in c(1:n))
  {
    m <- expansion(m, e)
    m <- inflation(m ,r)
  }

  print("End")
  print(m)

  return(m)
}

a <- matrix(c(0,1,1,0,1,0,0,0,0,0,0,1,0,0,1,0),4,4)
b <- mcl(a,2,2,10)

print(b)
